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Understanding stem cell (SC) population dynamics is essential for developing models that can be used in 
basic science and medicine, to aid in predicting cells fate. These models can be used as tools e.g. in studying 
patho-physiological events at the cellular and tissue level, predicting (mal)functions along the 
developmental course, and personalized regenerative medicine. Using time-lapsed imaging and statistical 
tools, we show that the dynamics of SC populations involve a heterogeneous structure consisting of multiple 
sub-population behaviors. Using non-Gaussian statistical approaches, we identify the co-existence of fast 
and slow dividing subpopulations, and quiescent cells, in stem cells from three species. The mathematical 
analysis also shows that, instead of developing independently, SCs exhibit a time-dependent fractal behavior 
as they interact with each other through molecular and tactile signals. These findings suggest that more 
sophisticated models of SC dynamics should view SC populations as a collective and avoid the simplifying 
homogeneity assumption by accounting for the presence of more than one dividing sub-population, and 
their multi-fractal characteristics. 

Stem cells are classically defined as unspecialized cells that can self-renew and give rise to differentiated cell 
types during embryogenesis, and in the adult, during tissue homeostasis or injury repair. These functions 
make them highly attractive to study for the purposes of understanding ontogeny and development, or for 
their potential use in regenerative medicine and tissue engineering. 

After more than 25 years of extensive research of numerous stem cell types, the field still struggles with how to 
define stem cells based on a molecular or chemical signature. Defining stem cells using molecular surface markers 
is a challenge. The lack of consistency in marker expression may be due the changing expression of markers 
during stem cell manipulation, or maturation, or to population heterogeneity. Technical differences between 
laboratories' methods and reagents can also contribute to challenges in defining stem cells based on markers. This 
study takes a system-level view on stem cells and particularly focuses on heterogeneity and population dynamics 
which are poorly understood and contribute to ambiguity in the identification of cells responsible for specific 
functions. 

The notion of a stem cell population which is comprised of a network of cells with interacting functions is rarely 
considered ex vivo. In vivo, it is well established that stem cells reside within a niche or microenvironment 
consisting of different cell types that provide physical and chemical supportive factors. However, the in vitro study 
of stem cells often does not consider a niche environment. Rather, attempts to study stem cells have predomi- 
nantly focused on the isolation of purified subsets of cells with specific markers or functions' '". Yet, several 
reports suggest that a population level exists for various stem cell types including hematopoietic stem cells 
(HSCs), mesenchymal stem cells (MSCs)" '''and muscle stem cells'''"^''. In support of this, several groups have 
shown that an individual cell from a stem cell population can re-establish the heterogeneous parent 
population^'"^^. 

The basic science challenges with population heterogeneity subsequently lead to issues related to their use in 
regenerative medicine, e.g., in ensuring cell potency or predicting ex vivo expansion or growth rates. Producing 
therapeutic doses of stem cells by ex vivo expansion requires what the FDA terms 'more-than-minimal manip- 
ulation'^'''^%hich carries with it the risks of stem cells becoming contaminated, genetically transformed, or 
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Figure 1 | Heterogeneity in Stem Cell Populations, (a) Stem cell populations are comprised of mitotically active (dividing) and mitotically inactive 
(nondividing: quiescent cells, terminally differentiated cells, senescent cells, and dying cells), (b) Proliferative heterogeneity results from asynchrony in life 
stage which leads, in part, to the presence of these subpopulations. (c) More specifically, asymmetric divisions and fates or symmetric divisions'"' lead to 
different subpopulations. Intrinsic or extrinsic factors play a role in whether the daughter cells are different due to internal cues or environmental cues, 
(d) An example of development of heterogeneity from intrinsic cell differences. Asymmetry in DNA strand segregation would allow for stem cell self- 
renewal. The (blue) stem cell would retain the oldest DNA strands. This cartoon shows segregation of 1 chromosome: the oldest/grandparent strand is 
blue and designated 1.0, and the parent strand (a copy of the grandparent strand) is red and designated 1.1. All other copies are dashed lines and 
designated copy numbers are 1.1.1, 1.1.1.1, 1.1.1.1.1 etc. If non-random strand segregation occurs among all chromatids in the cell, the result is 
asymmetric divisions and self- renewal of the stem cell. We show here how heterogeneity could result in the expanding population. If distinct phenotypes 
occur based on the DNA strand copy numbers, then these phenotypes can be categorized based on the template and copy number. E.g., after the 3 
divisions shown in the lineage tree above, there would be 1 stem cell (1.0 and 1.1 strands, pOphenotype), 3 pl cells would have (1.1/1.1.1 stands), 3 pll cells 
(1.1.1/1.1.1.1) and 1 pIII cell (1.1.1.1/1.1.1.1.1). 



functionally changed. Bio-manufacturing methods must predict the 
time required to obtain potent dose(s) of stem cells, yet minimize the 
amount of time that cells are manipulated ex vivo. Indeed, models 
which can accurately predict the growth rate of a heterogeneous 
population will be valuable tools in the development of a manufac- 
turing process that minimizes cell culture time and reduces exposure 
to foreign materials. Until now, very few approaches examine non- 
linear behavior of stem cell growth^" Rather, the basic exponential 
model which is used extensively in cell biology assumes a constant 
division time, and that all cells are dividing. As such, the proliferative 
heterogeneity of stem cell populations has only been addressed 
superficially by segregating populations into dividing and nondivid- 
ing cells in compartment models'" structure population mod- 
gjj5,32.42-49^ and agent-based models™'''. Few have addressed the 
potential existence of distinct dividing subpopulations within the 
heterogeneous stem cell population. For example, Glauche et al.'^ 
developed a nonlinear, adaptive model which accounts for two 
functional states-quiescence and proliferative-to explain HSC 
Bromodeoxy-Uridine (BrdU)-label dilution data. Wilson et al.'' 
and Foudi et al.""* show that similar label dilution data can be modeled 
by a two subpopulations with distinct division rates. 

Mathematical models must consider population heterogeneity 
that may manifest in individual cell proliferative behavior, in 
molecular activity and metabolism, and/or in cell morphology. In 
terms of proliferative heterogeneity, stem cells may progress through 
the cell cycle at different rates, or may exit the cell cycle to various 



fates (Fig. la). Stem cell populations consist of mitotically active 
(dividing) and mitotically inactive (nondividing) cells. Included in 
the non-dividing fraction are 1) quiescent cells, which have the theor- 
etical potential to re-generate the entire population, 2) differentiated 
cells, and 3) senescent cells. AU these subpopulations are present in a 
typical asynchronously growing population of cells (Fig. lb). The 
dividing fraction can be further distinguished into cells that give 
rise to symmetric and asymmetric divisions/fates' " (Fig. Ib&c). 
Asymmetric divisions/fates lead to population heterogeneity, which 
can be characterized using different parameters including molecular, 
morphological or proliferative features. Asymmetry may result from 
a number of mechanisms including the immortal strand hypothesis'' 
(Fig. Id). Proliferative heterogeneity, which we examine in this 
paper, results from the presence of all these subpopulations and 
asymmetric cell fates'". 

In this paper, we demonstrate that the combination of statistical 
tools and time-lapsed imaging of individual cells' dynamics allows us 
to study stem cells as a population or network of cells. This com- 
bination of imaging-based approaches as previously used by 
Bahnson et al."" or, more recently, Scherf et al." and Rappaport et 
al.'", with rigorous multivariate statistics, and data-driven modeling 
has been recognized recently as being of fundamental importance for 
the future of cell biology". As such, we hypothesize that higher level 
cell-cell interactions are at work to coordinate cell fate decisions 
leading to re-establishing heterogeneity. Indeed, we show here that 
multiscale phenomena do exist for individual cells in the population 
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Figure 2 | Time-lapsed image analysis of stem cell populations show heterogeneity in proliferative activity, (a) Cell division occurring in a human 
muscle stem cell population. Total time lapsed is 4.5 hours, scale bar represents 50 microns, (b) Individual cell tracking (Cytotracker,) through 4 
generations. The paths for parent and daughter cells are shown as an overlay on phase contrast images. Division times are used to construct cell lineage 
trees. The vertical distance on the corresponding lineage trees represents the cell division time or the amount of time the cell was on screen. 
Ipar; c) Heterogeneity in cell division times can be detected by cell lineage analysis (human muscle stem cells are shown here). Individual cells were 
examined using time-lapsed imaging over 3 days. Lineage history trees show the presence of both dividing and nondividing cells in the population. 
Variation in cell division times is also detected. Most current growth models assume homogeneity with a constant division time, and exponential growth, 
(d) Identification of dividing and nondividing cells in time-lapsed images after addition of BrdU to media for 48 hours and immunostaining for BrdU. 
Muscle stem cell populations contain both dividing (BrdU + cells, and nondividing (BrdU-) cells. Desmin staining illustrates that molecular 
heterogeneity for the myogenic cytoplasmic protein does not correspond with proliferative heterogeneity. 



for certain cell characteristics, namely ceU division time (DT); this 
process is power-law in nature and provides evidence for interactions 
between subpopulations. Consequently, in this study, we show that 
the dynamics of stem cell populations displays a heterogeneous 
structure consisting of multiple sub-population behaviors and sug- 
gest that rather than defining stem cells by exclusive cell markers, 
efforts should be made to understand the stem cell population, its 
subpopulations, and the interactions among them. 

Results 

Stem cell dynamics display a heterogeneous structure. Using time- 
lapsed imaging'"'"', we recorded the cell growth and examined the 



structure and cell division dynamics for both muscle stem ceUs'^ ''' 
from mouse and rat, and human mesenchymal stem cells" ''^. Time- 
lapsed images were acquired at 10-minute intervals over the course of 
5 days. From these images, we obtained direct measurements of cell 
cycle times (division times) as time elapsed between cytokinetic 
events (Fig. 2a). Based on these data, we were able to construct cell 
lineage trees of divisional history (Fig. 2b). 

Heterogeneity was detected in cell lineage analysis of all stem cells 
(Fig. 2c). Lineage history trees show the presence of both dividing 
and nondividing cells in the population. Variation in DTs is also 
detected. Proliferative heterogeneity was also identified by BrdU 
labeling as we observed both BrdU+ and BrdU- cells; this did not 
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correlate well with molecular heterogeneity in expression of the cyto- 
plasmic myogenic protein desmin (Fig. 2d), which illustrates the 
complexity of understanding heterogeneous characteristics. 

For the stem cell DTs from all populations examined, we first 
estimate the empirical cumulative distribution function (see 
Figure 3. a and Supplementary Material Figures l.a, 2. a, and 3. a), 
the empirical probability of observing a cell DT larger than a given 
threshold (also called exceedance probability), then perform a max- 
imum likelihood estimation (MLE) of well-known distributions and 
investigate their goodness-of-fit with respect to the actual data (see 
Supplementary Material Figures 1, 2, 3 and 4 and Notes 1, 2 and 3). 
The main reason for analyzing the exceedance probabilities is to 
quantify whether the stem cell growth exhibits an exponential or a 
power law type of behavior. By analyzing the exceedance probability 
of cell DTs for all three species, we find that the cell DTs are not well 
approximated by exponential type of distributions (Supplementary 
Material Figures, lb, 2b, 3b) as one would assume. This is also shown 
by the deviations of the empirical exceedance probabilities from the 
Gaussian (Figures 3.a-b), log-normal (Figure 3.c) and gamma 
(Figure 3.d) distributions. Instead, the graphical analysis in 
Figures 3.e-f, and the Supplementary Figures 1, 2 and 3 and the 
statistical results in Supplementary Tables 1, 2 and 3 show that the 
exceedance probabilities are better fitted by long-tail distributions 
such as cx-stable, Student-t and generalized extreme value (Gev) dis- 
tributions. However, to better assess the validity of modeling stem 
cell DTs via uni-modal distributions, we estimate the statistical good- 
ness-of-fit for each MLE fit which implies analyzing the fitting errors 
between the empirical and postulated probability density function. 
These results suggest that although the stem cell DTs exhibit a long- 
tail type of behavior, this cannot be well fitted by uni-modal distribu- 
tions alone (see Supplementary Notes 1, 2 and 3, and Tables 1, 2, 
and 3). 

A similar conclusion can be reached by analyzing the graphical 
fittings in Figures l.c, 2.c and 3.c in Supplementary Material of 
several well known unimodal distributions (i.e., Gaussian, Gamma, 
Log-normal, Generalized extreme value, Student-t, WeibuU, asym- 
metric a-stable distribution) for the mouse and human cell DTs. In 
addition to graphical evidence, the Kolmogorov-Smirnov test, the 
Akaike Information Criterion and Hartigan's Dip test'"'' (see Notes 4, 
5, and 6 in the Supplementary Material) suggest that stem cell DTs 
cannot be modeled by an uni-modal distribution. This suggests that 
there actually may exist two dividing subpopulations of stem cells, 
each characterized by distinct dynamics of stem cells division rate; 
this is reflected by different probability density functions (PDFs) of 
the DTs. 

To overcome the poor fitting of uni-modal distributions and 
explain the observed heterogeneity in the PDFs of stem cell DTs, 
we consider the bi-modal statistical investigation and two PDF can- 
didates: i) a bi-modal Gaussian distribution (Figure 4.a) and ii) a bi- 
modal asymmetric a-stable distributions (Figures 4.b-d). By 
employing the MLE strategy for a mixture of two asymmetric oc- 
stable distributions and selecting the best fit according to the best 
Kolmogorov-Smirnov (K-S) test result, we observe that the empirical 
PDFs of mouse, rat, and human stem cell DTs are well fitted by the bi- 
modal a-stable distribution. Indeed, both graphically and statistically 
by looking at the error bars, we can see that the fitting involving the 
a-stable distribution captures much better the overall trend than the 
Gaussian one (Figure 4.a). 

In sum, the magnitude of p -values of the K-S tests imply not only 
that we cannot reject the bi-modal a-stable distribution as a viable 
model for stem cell DTs, but also that the stem cells population 
exhibits a heterogeneous structure consisting of at least two sub- 
populations, namely a sub-population of stem cells that divide faster 
(with an average DT of approximately 12 hours and corresponding 
to first peak in the PDF plot of Figure 4.b), and a second sub- 
population of stem cells with a slower dynamics (with an average 



DT of approximately 18 hours corresponding to the second peak in 
the PDF plot of Figure 4.b) which exhibits a long-tail type of 
behavior. 

Stem cell division times exhibit non-stationary behavior. Besides 
the heterogeneous structure of stem cells population, we also observe 
that the empirical PDF estimated from stem cell DTs exhibits a 
pronounced time dependent behavior (Supplementary Material 
Figure 5). This time dependent behavior is illustrated not only by 
the changes in the shape of the empirical PDF at different time 
windows, but also by the estimated values of the bi-modal PDF. 
Indeed, the PDF investigation shows that stem cell growth is not a 
stationary process, but rather a highly dynamic one characterized by 
time-dependent PDF variation. 

To further investigate the existence of a non-stationary behavior, 
we measure the mean (Figure 5.a), variance (Figure 5.b), skewness 
(Figure 5.c), and kurtosis (Figure 5.d) over a sliding window of 80 cell 
DT recordings (from the time-ordered mouse stem cell DT data). 
Additionally, Figures 5.e-g show the error bars in the time variation 
of mean, variance, skewness, and kurtosis for mouse, rat, and human 
cells. We observe that these higher order moments vary with time 
which suggests that the cell division is a non-stationary process (see 
Supplementary Materials Figures 6 and 7). For instance, the mean of 
the rat stem cell DTs exhibit a wide range variation from 15 hours in 
the beginning to almost 19 hours at later times. Similarly, the vari- 
ance, skewness, and kurtosis exhibit a spiky behavior which cannot 
be modeled if one assumes a stationary behavior. In addition, the 
non-zero values for skewness and kurtosis show that the rat stem cell 
DTs do not display a Gaussian behavior; this supports our predic- 
tions concerning the existence of a non-Gaussian dynamics of stem 
cell growth shown in Fig. 3. Further, the non-zero skewness shows 
that the stem cell DTs are not symmetrically distributed around the 
mean and may exhibit different tails at both ends of the PDF curve. 
Similarly, the non-zero kurtosis implies that the stem cell DTs can 
exhibit large values and so a better modeling approach than the 
Gaussian-based framework needs to rely on long-tail distributions. 

Moreover, the observed variability and spiky dynamics of these 
higher moments imply that the stem cell DTs exhibit a heterosce- 
dastic dynamics (i.e., some sub-populations exhibit different mo- 
ment variability compared to others). An alternative proof of 
existence of heteroscedastic dynamics is shown in the Supplemen- 
tary Material Figure 8 where shuffled stem cell DT series lose their 
correlation structure when compared to the initial data sets. In gen- 
eral terms, the heteroscedastic dynamics means that the statistical 
moments of the stem cell DTs and, implicitly, those associated with 
stem cell growth exhibit a local variation which is very different from 
the global one or asymptotic limits. In addition, the very existence of 
this heteroscedastic behavior may not only confirm the heterogeneity 
of stem cell population, but also be indicative of an aging phenom- 
enon that is characteristic to all multi- cellular organisms. As we later 
show in the paper, this aspect can play a crucial role not only in 
constructing an accurate mathematical model of stem cells growth, 
but also for predicting the structure and size of the entire population 
over time with high confidence while accounting through non-sta- 
tionary methods and tools for the age of a patient. 

Stem cell growth rates possesses multi-fractal characteristics. For a 

comprehensive investigation of the heteroscedastic dynamics of stem 
cell growth, we investigate the relationship between the higher order 
moments of stem cells dynamics and their order; we also estimate 
both the multi-fractal spectrum and generalized Hurst exponent 
function [A fractal is a geometrical object or stochastic process that 
displays self- similarity, on all scales and is characterized by a single 
fractal dimension. The fractal dimension is a mathematical concept 
measuring the unique features such as geometrical shape of an object 
or irregularity of a stochastic process. The object need not exhibit 
exactly the same structure at all scales, but the same "type" of 
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Figure 3 | Statistical investigation and interpretation of tracked stem cell division times, (a) Comparison between the cumulative distribution function 
(CDF) and the 95% confidence intervals for the Gaussian distribution fitted via maximum hkelihood estimation method and the empirical CDF of the 
mouse stem cell division times, (b) Comparison between the exceedance probability (i.e., probability of observing a stem cell division time greater than a 
specific threshold) and the 95% confidence interval for the fitted Gaussian distribution and the empirical inverse CDF of the mouse stem cell division 
times, (c) Comparison between the exceedance probability and its 95% confidence intervals for the fitted Log-normal distribution and the empirical 
inverse CDF of the mouse stem cell division times, (d) Comparison between the exceedance probability and its 95% confidence intervals for the fitted 
Gamma distribution and the empirical inverse CDF of the mouse stem cell division times, (e) Maximum likelihood fitting and the error bars for 
the uni-modal Gaussian probability density function (PDF). The empirical PDF is shown by the blue dots, (f) Maximum likelihood fitting and the 
error bars for the a-stable PDF. Rejection of uni-modal mathematical modeling approaches for stem cell growth is motivated not only by the 
inconsistencies observed via graphical inspection of the empirical PDF against postulated PDF, but also by the very small p-value probabilities of the 
Kolmogorov-Smirnov (K-S) test (see the Supplementary Material). 
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Figure 4 | Statistical analysis of bi-modal asymmetric a-stable distribution fitting and multi- fractal analysis of stem ceU division tim (a) Comparison 
between the empirical PDF (black dots) of mouse stem cell division times (DTs) and the bi-modal Gaussian distribution (red line), (b) Comparison 
between the empirical PDF of mouse stem cell DTs and the bi-modal a-stable distribution for all 889 samples. The high p-value of 0.78 for the K-S test 
shows that the postulated bi-modal ot-stable distribution cannot be rejected as a model, (c) Comparison between the empirical PDF of rat stem cell DTs 
and the bi-modal a-stable distribution for all 1055 samples. The p-value of 0. 11 for the K-S test shows that the postulated bi-modal a-stable distribution 
cannot be rejected as a model, (d) Comparison between the empirical PDF of human MSG DTs and the bi-modal a-stable distribution for aU 350 samples. 
The p-value of 0.8 1 for the K-S test shows that the postulated bi-modal a-stable distribution cannot be rejected as a model, (e) Multi-fractal spectrum as a 
function of Lipschitz-Holder mass exponent for mouse, rat and human muscle derived stem cell division times, (f) Generalized Hurst exponent as a 
function of q-th order moment for mouse, rat and human muscle derived stem cell division times. 



structures must appear on all scales. Similarly, a stochastic process is 
regarded as fractal if its probability distribution function measured at 
two scales is self-similar and related to each other via a power law 
relationship. In contrast, a multi-fractal object or stochastic process 



is characterized by a series of fractal dimensions.]. The rationale for 
estimating the generalized Hurst function in Figure 4.f is two-fold: 
( 1 ) It allows us to detect the existence of multi-fractality. For instance, 
if the generalized exponent is independent of the order q of the 
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Figure 5 | Higher order moment analysis of mouse, rat, and human stem cells division times, (a) Moving average analysis of mean of mouse stem 
cell division times over a time window of 80 samples, (b) Moving average analysis of variance of mouse stem cell division times over a time window of 80 
samples, (c) Moving average analysis of skewness of mouse stem cell division times over a time window of 80 samples, (d) Moving average analysis of 
kurtosis of mouse stem cell division times over a time window of 80 samples. Moving average analysis of higher order moments of stem cell division 
times shows that cell growth exhibits a pronounced time dependent behavior, (e) Average values and errors for the mean, variance, skewness and kurtosis 
of the mouse stem cell division times, (f) Average values and errors for the mean, variance, skewness and kurtosis of the rat stem cell division times, 
(g) Average values and errors for the mean, variance, skewness and kurtosis of the human MSG division times. 



moment (see Figure 4.f), then the stochastic process is mono-fractal. 
In contrast, if the generalized Hurst exponent is dependent of the 
order q of the moment, then the process is called multi-fractal (see 
Figure 4.f). (2) It allows us to compute the Hurst parameter and 
discriminate between short-range and long-range memory effects. 
For instance, if the Hurst exponent equals 0.5 then we can state that 
the dynamics of stem cell populations is governed by short-range 
memory or Markovian models. In contrast, if the Hurst exponent 
ranges between 0.5 and 1, then the dynamics of stem cell populations 
is said to possess long-range memory characteristics and non- 
Markovian dynamical models are needed. While mono-fractality 
implies that a single Lipschitz-Holder exponent*"' (i.e., a function 
f{t) is said to have Lipschitz-Holder exponent a around point t if 
f(t + x) — f(t) ~ x") characterizes the entire cell dynamics, the multi- 
fractal behavior shows that the existing temporal heterogeneity in 
stem cell population leads to multiple Lipschitz-Holder exponents. 
The multi-fractal spectrum encapsulates in its distribution of 
Lipschitz-Holder exponents information about multiple point 
correlations existing both in space and time. Implicitly, it is also a 
measure of the memory characterizing the stem cell dynamics. 

By measuring the dependency of the higher order statistical 
moments on their order value and using the Legendre transforma- 
tion, we estimate the multi-fractal spectrum for the mouse, rat, and 
human stem cell DTs (Figure 4.e and Supplementary Material Figure 
9). One can clearly observe from Figure 4.e not only that the stem cell 
DTs exhibit a wide multi-fractal spectrum, but also the Lipschitz- 
Holder exponents are mostly negative. The fact that the stem cell DTs 
display a significant range of negative Lipschitz-Holder exponents 
implies that the stem cell growth process exhibits an oscillating 
behavior consisting of active periods of growth (when the majority 
of the cells are actively dividing) followed by periods when fewer cells 
are recruited into the cell cycle. Alternatively, the results in Figure 4.f 
(and Supplementary Material Figure 10) display a complex nonlinear 



dependency between the generalized Hurst exponent and the q-th 
order moment, which again supports the existence of multi-fractal 
behavior in stem cell growth. 

Based on these findings, our subsequent question is whether the 
current models relying on the assumption of either having constant 
or exponential growth in stem cell population size are suitable for 
predicting the population size; answering this question is critical for 
designing cell therapies or tissue engineered constructs for future 
personalized medicine. Therefore, we next considered the time 
lapsed traces of the number of stem cells in the population over time 
and investigate the goodness-of-fit for three mathematical modeling 
approaches: /) exponential growth associated with a non-fractal 
dynamics (Figure 6.a); »') power law growth corresponding to a 
mono-fractal dynamics (Figure 6.b); and Hi) a stretched exponential 
hinting towards a multi-fractal approach (Figure 6.c). We observe 
that the stretched exponential and power law type of models offer 
better prediction capabilities than the state-of-the art exponential 
models (see Figure 6.d). By ignoring the multiscale biological 
dynamics and spatio-temporal correlations, the exponential models 
can be misleading. Indeed, the exponential growth assumption can 
lead to very optimistic predictions concerning the population size of 
stem cells and this may have detrimental effects on the therapeutic 
strategies. Our results suggest that the multi-fractal behavior of stem 
cell dynamics can be modeled as a superposition of several stretched 
exponential functions. 

Discussion 

Here, we show for the first time that distinct subpopulations of 
dividing cells exist within several heterogeneous stem cell popula- 
tions. We report that the populations show time-dependent and 
fractal behavior. Both newborn derived mesenchymal stem cells 
and adult derived muscle stem cells demonstrate these population 
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Figure 6 | Implications of multi-fractal behavior on population size and their contrasting differences when considering simplified models. 

(a) Comparison between the empirical average number of rat stem cells as a function of time and the exponential growth model, (b) Comparison 
between the empirical average number of rat stem cells as a function of time and the power-law growth model, (c) Comparison between the empirical 
average number of rat stem cells as a function of time and the stretched exponential growth model, (d) Empirical average number of rat stem cells as a 
function of time and three mathematical interpolations, i.e. exponential (red line), power law (black dotted line), and stretched exponential (green line). 
Their best parameters are shown in the legends in plots a-c. The first 75% of the data samples ( ~ 150 hours) also shown in black circles are used to identify 
the parameters of each model. The remaining 25% of the data samples (—50 hours) are used to test the prediction capabilities of each model. One can 
notice that the stretched exponential and power law models have better prediction capabilities than the simple exponential model. This is in agreement 
with our multi-fractal analysis which predicts that the intrinsic growth of stem cell population is characterized by non-exponential features. 



characteristics. By employing non-Gaussian statistics, statistical 
hypothesis testing and statistical inference, we identify the presence 
of two unique dividing subpopulations of stem ceils, each with a 
distinct cell division rate and characterized by a different probability 
density function of DTs. Previous reports treated stem cells as homo- 
geneous populations and assumed either a constant division time for 
the whole population or Gaussian distributed DTs. In terms of 
improving how stem cells are defined, our study strongly supports 
the notion of defining cells from a population level perspective in 
order to account for cell-cell interactions that affect the population 
structure (composition of different subpopulations) and dynamics 
(interactions of and changes in the subpopulations). 

Bacteria and yeast cell division and growth rates have been exam- 
ined extensively for more than 50 years; the distributions appear to be 
rather homogeneous'"'''" and, to the best of our knowledge, there are 
no reports of bi- or multi-modal growth rates. For instance, Tsuru et al.™ 
recognize the existence of fluctuations in E.coli growth rates and 
investigated how the changes in cell volume of individual bacteria 
are correlated with fluctuations in protein concentration. Comple- 
menting such a phenotypic diversity, a report by Niven et al. descri- 
bes skewed cell division distributions in £. coM in the presence of 
environmental factors''. The authors examined several distributions 
to fit their data (also obtained via time-lapsed cell imaging) and no 



definitive model distribution was identified, possibly due to the small 
sample size''. 

More recent studies of differentiated cells, particularly those that 
use similar technology to permit the collection of a dataset of indi- 
vidual cell measurements, also did not report bi- modality"''"''''''''''^. For 
instance, Tyson et al." emphasized the importance of investigating 
the statistics of intermitotic times and modeled these times via an 
exponentially modified Gaussian distribution. Based on their ana- 
lysis, Tyson et al. concluded that their proposed exponentially 
modified Gaussian models are not definitive and require further 
improvement to accurately accommodate age structure and 
dynamics of stem cell populations. Along the same lines, Chang 
et al.'' demonstrate the existence of cell-to-cell variability in clonal 
populations of mouse hematopoietic progenitor cells and showed a 
link to heterogeneity of gene expression. Tan et al.'" used single-cell 
gene expression to show functional heterogeneity in epidermal cells, 
while Koschmieder et al.™ and Cantor et al."" investigated it in the 
hematopoietic lineage. We see that such mechanism might be an 
underlying cause of the heterogeneity in our populations of stem 
cells. Furthermore, in our own previous reports, we used both stand- 
ard Gaussian and non-parametric statistics to examine the distribu- 
tions of cell cycle times in human and mouse stem cells. However, 
lack of sophistication in our previous statistical analyses and the 
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small size^''", did not cillow us to identify bi-modality (let alone 
multi-modality) in stem cells dynamics which we report here. 

Stem cell population heterogeneity may develop through the pro- 
cess of spontaneous progression of stem cell differentiation. It is 
worth mentioning that the interaction of the cells with the culture 
substrate and the rigidity of the substrate are important factors in 
modulating differentiation and cell phenotype heterogeneity as has 
been shown by Gilbert et al."'. In most scenarios, non-dividing qui- 
escent cells are presumed to become activated to slowly dividing 
multipotent stem cells which then progress to transiently amplifying 
(TA) committed progenitors and finally non-dividing differentiated 
cells. Here, we detected, in all populations the presence of a non- 
dividing quiescent subpopulation that could re-enter the cell cycle. 
The slowly dividing subpopulation for MSCs (16.8 hrs), rat stem 
cells (21.2 hrs), or mouse stem cells (21.9 hrs) was at least 4 hours 
slower than the fast diving subpopulation (human MSG 12.8 hrs, rat 
stem cells 11.7 hrs, mouse stem cells 13.2 hrs). It is possible that the 
slow dividing fraction represents multipotent cells, while the faster 
dividing cells represent committed precursors or transiently-amp- 
lifying cells. However, at this point, there are no clear methods to 
physically separate these subpopulations, and BrdU pulsing experi- 
ments (Fig. 2) did not show a correlation between the myogenic 
precursor marker Desmin and BrdU rate of incorporation in stem 
cell populations. Nevertheless, the statistical analysis demonstrates 
the presence of two distinct populations which grow at different 
rates. This was detected in different species, and in different 
stem cell types, which suggests that it may be indeed a universal 
characteristic. 

The existence of these different phenotypes may allow for both 
short-term and long-term functions of the stem cell populations in 
response to tissue injury. Indeed, fast dividing TA cells can provide 
immediate restoration of damaged/lost differentiated cells, while the 
dynamic of TA cell depletion for acute repair may illicit a faster 
transition rate for cells in the multi-potent state to move to the TA 
state to re-establish some equilibrium of subpopulations. In this 
scenario, because dividing cells expand according to a stretched 
exponential, the quiescent stem cell activation could be conservative 
and occur only in the case of severe tissue loss or damage. In the case 
of normal tissue homeostasis, we could expect to observe oscillating 
dynamics (Figure 5. a or Supplementary Material Figures 6.a-c and 
7.a-c) related to cell division behavior in the population, as normal 
cell death occurs, and lost cells are replaced. 

Further, our statistical analysis demonstrates that the population 
subsets display a non-stationary (i.e., time-dependent statistics). 
Non-stationary statistical analysis shows that not only do individual 
cells change over time, intrinsically, but the stem cell population also 
changes over time. In part, the non-stationary phenomenon finding 
shows that stem cells exhibit an aging process which is characteristic 
to all living organisms. Gonsequently, the rates of transition between 
subpopulation states will be important for understanding their 
evolution over time. 

The fractal behavior we observe in stem cell DTs identifies a 
stochastic process that displays a probability density function that 
is self-similar in nature**^ *', i.e., the distribution of cell DTs at one 
scale can be retrieved from that of another scale by using the fractal 
dimension and a scaling coefficient. Simply speaking, the existence of 
fractal statistics in stem cell DTs implies that the birth of a new cell is 
not a random (i.e., uncorrelated) event from the development of 
previous ones. In other words, there is some form of correlation 
between cell divisions and this leads to a complex behavior which 
cannot be quantified by current mathematical models of stem cell 
growth. Nevertheless, due to numerous sources of variation (e.g., 
intrinsic differences, environmental interactions, nutrient availabil- 
ity, cell density), stem cell DTs cannot be characterized by a single 
fractal dimension, but by a set of scaling exponents which also 
demonstrates the existence of heterogeneity in stem cell population 



structure. Taking into account such characteristics allows us to con- 
struct an accurate dynamical system approach as described by 
Furusawa and Kaneko"'''''^ which may explain and, more importantly, 
predict the stem cell dynamics. Furthermore, the newly proposed 
model can be extended to include cell cycle and cell volume, along 
the lines of Halter et al.""" and Anderson et al."'; these parameters may 
be used to show the multi-fractal behavior as the cells interact with 
each other. 

The value for teasing out the subsets of cells, their interactions, 
time-dependent and fractal behavior has implications for applied cell 
therapeutics. First, it is needed to identify/define the stem cell prod- 
uct. Second, it is necessary to accurately predict ex vivo growth rates 
of the stem cell populations. In regards to developing an FDA- 
approved cell therapeutic, purity and potency criteria need to be 
established. Recognizing that subpopulations, and hence multiple 
targets, exist in the stem cell population is essential to developing 
tools to control the fate of stem cells. Although, we and others have 
shown that the quiescent population exists and requires nonlinear 
growth models^" our study is the first to show quantitatively and 
qualitatively that more sophisticated growth models are warranted 
due to the presence of multiple phenotypes. Models which account 
for population dynamics can better predict whether cytokines or 
growth factors that increase proliferation rates (e.g., for the purposes 
of obtaining certain cell doses for transplantations) cause perturba- 
tions in the dynamics of subpopulations, and how this affects the 
behavior of the entire cells population. 

This statistical approach may also be exploited for personalized 
medicine strategies. Using high throughput live cell imaging and our 
computational analyses, inter-patient variability and intrinsic intra- 
patient variation in stem cells dynamics can be studied. In the era of 
individualized medicine, the computational modeling may also 
include genomics and proteomics information and the dynamics 
of an individual biopsy-derived stem cell growth can be assessed 
before developing the cell therapy approach. Personalized medicine 
could also exploit our statistical investigation for estimating the rich- 
ness of fractal behavior at the cellular level and make predictions at 
higher scales. By early detection of losing fractal richness at the level 
of stem cell evolution, one could predict the signs of genetic instab- 
ility and onset of a disease well in advance of actual symptoms in the 
patient, when it may be hard, or impossible, to intervene. 

Further studies are needed to examine the molecular phenotype 
and transition rates of the two dividing subsets identified here. 
Epigenetic modifications may account for differences in phenotype 
or division rates, or the ceU-cell communications may modulate 
division rates. The oscillating time-dependent behavior together with 
the finding that these subpopulations interact with each other (i.e., 
cells do not grow in isolation, but instead interact at various scales). 
In addition to that, the cells indeed interact with their microenviron- 
ment^' including the extracellular matrix"', growth factors, chemo- 
kines'", and the substrate they are adhering to*'. These interactions 
have been shown to affect the differentiation and fate decisions of 
stem cells. All of these factors call for developing new mathematical 
models describing their interaction and overall population growth. A 
challenge to these studies will be the isolation of these subsets and 
maintenance without re-establishing the heterogeneous parent 
population. 

Methods 

We used time-lapsed imaging*'"'^^ to directly examine the cell growth of mouse and rat 
stem cells isolated from skeletal muscle tissue and human mesenchymal stem cells. 
(Additional detail regarding imaging experimental set-up is included in''" ''^ and stem 
cell isolation'''^ ''^ is included in the Supplementary Material Note) .Time-lapsed images 
were acquired at 10-minute intervals over the course of 5 days on collagen-coated 
surfaces. To determine the divisional status of the cells, and the cell cycle length of 
dividing cells, we used a custom developed program, i.e., ImageView {Karios 
Instruments, Harmar, PA). Cell cycle length was recorded as the time elapsed between 
cytokinetic events, as previously reported^^''"''^'. Dividing cells were classified as cells 
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with a visible DT. Non-dividing cells were classified as cells that were visible on screen 
for over twice the average DT. 

After determining DTs, the cells were tracked in order to perform lineage analysis 
and determine differences between dividing and non-dividing cells while also iden- 
tifying the quiescent cells. A custom built image analysis program, CytoTracker, 
(Karios Instruments, Harmar,PA), was used to perform lineage analysis. A cell was 
selected and followed for the entire time it was visible on screen. When a cell divides, 
the division event is recorded, and each daughter cell was subsequently tracked. The 
entire cell lineage is tracked until either the end of the video or the cells are lost to 
follow up. Each cell is assigned a name in order to keep track of the lineage. 

The initial parent cell is 1.0. The two daughter cells are named 1.1 and 1.2. Their 
daughter cells are named 1.11,1.12 and 1.21, 1.22 respectively. The tracks/locations {x 
and y location) for the entire family of cells are recorded along with the cell name for 
each jpeg image.This allows a cells lineage to be determined. The cell histories can 
then be converted to a graphical display, or lineage tree which displays all the cells in 
the family, and their DTs, or their time on screen if the cell was lost to follow up 
(Figure 2). In the graphical displays, the vertical lines represent a cell, and the length of 
the line corresponds to how long the cell was on screen. Lineage trees allow the 
divisional status of a cell to be determined. A tree with many cells represents a cell that 
is actively dividing and giving rise to progeny that continue to divide and contribute to 
renewal of the population. The lineage trees also allow DTs to be compared, or to 
identify cells that are non-dividing and visible on screen for over twice the division 
time. Also, asymmetrical divisions can easily be identified. 
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